/***************************************************************************
specfunc.c - SPECIAL FUNCTIONS FOR DSP


This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

Copyright (C) 2013  Bakhurin Sergey
***************************************************************************/

#define _USE_MATH_DEFINES
#include <math.h>

#include "specfunc.h"
#include "dsp.h"
#include "cmplx.h"




/***************************************************************************
INVERSE HYPERBOLIC SINE asinh(x) = log(x +sqrt(x^2+1.0))

Input parameters:
	typedsp x - input argument 

Return: 
	inverse hyperbolic sine of x */

typedsp asinh(typedsp x)
{
	return log(x +sqrt(x*x+1.0));
}




/***************************************************************************
HYPERBOLIC COSINE OF REAL ARGUMENT 

Input parameters:
	typedsp x - input argument 

Return: 
	hyperbolic cosine of x */

typedsp cosh(typedsp x)
{
	return 0.5*(exp(x) + exp(-x));
}




/***************************************************************************
 JACOBI ELLIPTIC cd FUNCTION

Input parameters:
	typedsp u - elliptic amplitude 
	typedsp k - elliptic modulus 
Return: 
	cd(u,k)*/

typedsp ellipCd(typedsp u, typedsp k)
{
	long n;
	typedsp ktmp[ELLIP_ITER_COUNT], w;
	ktmp[0] = k;
	for(n = 1; n <  ELLIP_ITER_COUNT; n++)
	{
		ktmp[n]  = ktmp[n-1]/(1.0+sqrt(1.0-ktmp[n-1]*ktmp[n-1]));
		ktmp[n] *= ktmp[n];
	}
	w = cos(u*M_PI_2);
	for(n = ELLIP_ITER_COUNT-1; n > 0; n--)
	{
		w = (1.0 + ktmp[n])/(1.0/w + ktmp[n] * w);
	}

	return w;
}




/***************************************************************************
 COMPLETE ELLIPTIC INTEGRL FUNCTION

Input parameters:
	typedsp k - elliptic modulus 
Return: 
	elliptic integral of argument k*/

typedsp ellipInt(typedsp k)
{
	long n;
	typedsp ktmp, K;
	ktmp = k;
	K = M_PI_2;
	for(n = 1; n <  ELLIP_ITER_COUNT; n++)
	{
		ktmp  = ktmp/(1.0+sqrt(1.0-ktmp*ktmp));
		ktmp *= ktmp;
		K *= (1.0+ktmp);
	}
	return K;
}




/***************************************************************************
 JACOBI ELLIPTIC sn FUNCTION

Input parameters:
	typedsp u - elliptic amplitude 
	typedsp k - elliptic modulus 
Return: 
	sn(u,k)*/

typedsp ellipSn(typedsp u, typedsp k)
{
	long n;
	typedsp ktmp[ELLIP_ITER_COUNT], w;
	ktmp[0] = k;
	for(n = 1; n <  ELLIP_ITER_COUNT; n++)
	{
		ktmp[n]  = ktmp[n-1]/(1.0+sqrt(1.0-ktmp[n-1]*ktmp[n-1]));
		ktmp[n] *= ktmp[n];
	}
	w = sin(u*M_PI_2);
	for(n = ELLIP_ITER_COUNT-1; n > 0; n--)
	{
		w = (1.0 + ktmp[n])/(1.0/w + ktmp[n] * w);
	}

	return w;
}




/***************************************************************************
POLYNOM VALUE CALCULATION FOR COMPLEX ARGUMENT x

P(x) = a[0] + a[1]*x + a[2]*x^2 + ... + a[n-1] x^(n-1)

Input parameters:
	typedsp *a - Pointer to polynom rel vector coeff  
	long n	   - "a" coeff count. Polynom max power is n-1
	complex x  - complex argument	
Return: 
	polynom vlue of complex argument "x"*/

complex polyValCmplx(typedsp *a, long n, complex x)
{
	complex  p;
	long k;

	p = numCmplx(a[n-1], 0.0);
	for(k = n-2; k > -1; k--)
	{
		p = sumCmplx(numCmplx(a[k], 0.0), mulCmplx(x, p));
	}
	return p;
}



/***************************************************************************
HYPERBOLIC SINE OF REAL ARGUMENT 

Input parameters:
	typedsp x - input argument 

Return: 
	hyperbolic sine of x */
typedsp sinh(typedsp x)
{
	return 0.5*(exp(x) - exp(-x));
}